#include <list>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <random>
#include <cassert>

#include "GeoGrid.h"

void GeoDelaunay::SaveToVTK(std::string filename) {
        // Open the VTK file for writing
    std::ofstream vtkFile(filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }
    printf("Delaunay Save VTK to %s\n", filename.c_str());

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET UNSTRUCTURED_GRID\n";
    vtkFile << "\n";

    int n_vertices = d->nb_vertices();
    const double* vertices = d->vertices_ptr();
    // Write the grid points to the VTK file
    vtkFile << "POINTS " << n_vertices << " float\n";
    for(int i = 0; i < n_vertices; i++) {
        vtkFile << vertices[3*i] << " " << vertices[3*i+1] << " " << vertices[3*i+2] << "\n";
    }
    vtkFile << "\n";

    int n_cells = d->nb_cells();
    int valid_cells = 0, invalid_cells = 0;

    for(int i = 0; i < n_cells; i++) {
        const GEO::signed_index_t* cell = d->cell_to_v() + 4*i;
        int flag = 0;
        for(int j = 0; j < 4; j++) {
            if(cell[j] < 0) {
                invalid_cells++;
                flag = 1;
            }
        }
        if(!flag) {
            valid_cells++;
        }
    }

    printf("valid Delaunay cells: %d, invalid Delaunay cells: %d\n", valid_cells, invalid_cells);

    // Write the grid cells to the VTK file
    vtkFile << "CELLS " << valid_cells << " " << valid_cells * 5 << "\n";
    for(int i = 0; i < n_cells; i++) {
        const GEO::signed_index_t* cell = d->cell_to_v() + 4*i;
        int flag = 0;
        for(int j = 0; j < 4; j++) {
            if(cell[j] < 0) {
                flag = 1;
            }
        }
        if(!flag) {
            vtkFile << "4 ";
            for(int j = 0; j < 4; j++) {
                vtkFile << cell[j] << " ";
            }
            vtkFile << "\n";
        }
    }

    // Write the cell types to the VTK file
    vtkFile << "CELL_TYPES " << valid_cells << "\n"; 
    for (int i = 0; i < valid_cells; i++) {
        vtkFile << "10\n";
    }

    // Close the VTK file
    vtkFile.close();
}

void GeoDelaunay::get_cell(GEO::index_t v, GEO::ConvexCell& C) {
    d->copy_Laguerre_cell_from_Delaunay(v, C, W_);
    if(!periodic) {
        C.clip_by_plane(GEO::vec4( 1.0, 0.0, 0.0, 0.0));
        C.clip_by_plane(GEO::vec4(-1.0, 0.0, 0.0, 1.0));
        C.clip_by_plane(GEO::vec4( 0.0, 1.0, 0.0, 0.0));
        C.clip_by_plane(GEO::vec4( 0.0,-1.0, 0.0, 1.0));
        C.clip_by_plane(GEO::vec4( 0.0, 0.0, 1.0, 0.0));
        C.clip_by_plane(GEO::vec4( 0.0, 0.0,-1.0, 1.0));
    }
    C.compute_geometry();
}

void GeoDelaunay::Lloyd_iteration() {
    std::vector<double> new_points(flat_points.size());
    GEO::ConvexCell C;
    for(GEO::index_t v=0; v < flat_points.size()/3; ++v) {
        get_cell(v, C);
        GEO::vec3 g = C.barycenter();
        new_points[3*v]   = g.x;
        new_points[3*v+1] = g.y;
        new_points[3*v+2] = g.z;
    }
    // In periodic mode, points may escape out of
    // the domain. Relocate them in [0,1]^3
    for(GEO::index_t i=0; i<new_points.size(); ++i) {
        if(new_points[i] < 0.0) {
            new_points[i] += 1.0;
        }
        if(new_points[i] > 1.0) {
            new_points[i] -= 1.0;
        }
    }
    flat_points.swap(new_points);
    d->set_vertices(flat_points.size()/3, flat_points.data());
    d->compute();

}

void GeoDelaunay2d::SaveToVTK(std::string filename) {
        // Open the VTK file for writing
    std::ofstream vtkFile(filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }
    printf("Delaunay2d Save VTK to %s\n", filename.c_str());

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET POLYDATA\n";
    vtkFile << "\n";

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << n_points << " float\n";
    for(int i = 0; i < n_points; i++) {
        vtkFile << flat_points[2*i] << " " << flat_points[2*i+1] << " " << 0 << "\n";
    }
    vtkFile << "\n";

    int n_trangles = d->nb_cells();

    printf("Delaunay trangles: %d\n", n_trangles);

    // Write the grid cells to the VTK file
    vtkFile << "LINES " << n_trangles * 3 << " " << n_trangles * 3 * 3 << "\n";
    for(int i = 0; i < n_trangles; i++) {
        const GEO::signed_index_t* cell = d->cell_to_v() + 3*i;
        for(int j = 0; j < 3; j++) {
            vtkFile << 2 << " " << cell[j] << " " << cell[(j+1)%3] << "\n";
        }
    }

    // Close the VTK file
    vtkFile.close();
}

void LayerSurface::SaveToVTK(std::string filename) {
        // Open the VTK file for writing
    std::ofstream vtkFile(filename);
    if (!vtkFile.is_open()) {
        // If the file cannot be opened, print an error message and return
        std::cerr << "Unable to open file: " << filename << std::endl;
        return;
    }
    printf("Delaunay2d Save VTK to %s\n", filename.c_str());

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET POLYDATA\n";
    vtkFile << "\n";

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << pts.size() << " float\n";
    for(int i = 0; i < pts.size(); i++) {
        vtkFile << pts[i].x << " " << pts[i].y << " " << pts[i].z << "\n";
    }
    vtkFile << "\n";

    int n_trangles = tris.size();

    // Write the grid cells to the VTK file
    // vtkFile << "LINES " << n_trangles * 3 << " " << n_trangles * 3 * 3 << "\n";
    // for(int i = 0; i < n_trangles; i++) {
    //     for(int j = 0; j < 3; j++) {
    //         vtkFile << 2 << " " << tris[i].v[j] << " " << tris[i].v[(j+1)%3] << "\n";
    //     }
    // }
    vtkFile << "POLYGONS " << n_trangles << " " << n_trangles * 4 << "\n";
    for(int i = 0; i < n_trangles; i++) {
        vtkFile << 3 << " " << tris[i].v[0] << " " << tris[i].v[1] << " " << tris[i].v[2] << "\n";
    }

    // Close the VTK file
    vtkFile.close();
}

double LayerSurface::TriPlaneIntersection(Triangle tri, Point3d p) {
    Point3d p0 = pts[tri.v[0]];
    Point3d p1 = pts[tri.v[1]];
    Point3d p2 = pts[tri.v[2]];
    Point3d n = geo3d::crossProduct(Point3d(p1.x - p0.x, p1.y - p0.y, p1.z - p0.z), Point3d(p2.x - p0.x, p2.y - p0.y, p2.z - p0.z));
    double t = p0.z - (n.x * (p.x - p0.x) + n.y * (p.y - p0.y)) / n.z;
    return t;
}

void LayerSurface::UpdateTriRange() {
    for(int i = 0; i < tris.size(); i++) {
        tris[i].min_brd = Point3d(1e9, 1e9, 1e9);
        tris[i].max_brd = Point3d(-1e9, -1e9, -1e9);
        for(int j = 0; j < 3; j++) {
            tris[i].min_brd.x = std::min(tris[i].min_brd.x, pts[tris[i].v[j]].x);
            tris[i].min_brd.y = std::min(tris[i].min_brd.y, pts[tris[i].v[j]].y);
            tris[i].min_brd.z = std::min(tris[i].min_brd.z, pts[tris[i].v[j]].z);
            tris[i].max_brd.x = std::max(tris[i].max_brd.x, pts[tris[i].v[j]].x);
            tris[i].max_brd.y = std::max(tris[i].max_brd.y, pts[tris[i].v[j]].y);
            tris[i].max_brd.z = std::max(tris[i].max_brd.z, pts[tris[i].v[j]].z);
            tris[i].coord_2d[j] = Point2d(pts[tris[i].v[j]].x, pts[tris[i].v[j]].y);
            tris[i].p[j] = pts[tris[i].v[j]];
        }
    }
}


int LayerSurface::checkInsertWithUpLine(Point3d p) {
    geo2d::Point2d p2d = geo2d::Point2d(p.x, p.y);
    Line up_line = Line(p, Point3d(p.x, p.y, p.z + 1e3));
    int res = 0;
    std::vector<Triangle*> check_tris = kdtree.IntersectLine(up_line);
    for(int i = 0; i < check_tris.size(); i++) {
        if(geo2d::isPointInsideTriangle(check_tris[i]->coord_2d[0], check_tris[i]->coord_2d[1], check_tris[i]->coord_2d[2], p2d)) {
            double intersect = TriPlaneIntersection(*check_tris[i], p);
            if(intersect > p.z) {
                res++;
            }
        }
    }
    return res;
}

int LayerSurface::checkInsertWithDownLine(Point3d p) {
    geo2d::Point2d p2d = geo2d::Point2d(p.x, p.y);
    Line up_line = Line(p, Point3d(p.x, p.y, p.z + 1e3));
    int res = 0;
    std::vector<Triangle*> check_tris = kdtree.IntersectLine(up_line);
    for(int i = 0; i < check_tris.size(); i++) {
        if(geo2d::isPointInsideTriangle(check_tris[i]->coord_2d[0], check_tris[i]->coord_2d[1], check_tris[i]->coord_2d[2], p2d)) {
            double intersect = TriPlaneIntersection(*check_tris[i], p);
            if(intersect < p.z) {
                res++;
            }
        }
    }
    return res;
}

Point3d ExpandInfo::GetTetrahedronPosition(const Point3d &a, const Point3d &b, const Point3d &c, const Point3d &d) {
    Point3d center = geo3d::computeMasscenter(a, b, c, d);
    Point3d ret;
    if(center.x < inner_rmin.x) {
        ret.x = -1;
    }
    else if(center.x > inner_rmax.x) {
        ret.x = 1;
    }
    else {
        ret.x = 0;
    }
    if(center.y < inner_rmin.y) {
        ret.y = -1;
    }
    else if(center.y > inner_rmax.y) {
        ret.y = 1;
    }
    else {
        ret.y = 0;
    }
    if(center.z < inner_rmin.z) {
        ret.z = -1;
    }
    else if(center.z > inner_rmax.z) {
        ret.z = 1;
    }
    else {
        ret.z = 0;
    }
    return ret;
}

GeoGrid::GeoGrid(GRDECLData* data) {
    GEO::initialize();
    GEO::CmdLine::import_arg_group("standard");
    GEO::CmdLine::import_arg_group("pre");
    GEO::CmdLine::import_arg_group("algo");

    this->data = data;
    
    printf("nx : %d, ny : %d, nz : %d\n", data->nx, data->ny, data->nz);

    // 添加非重复点
    ind_layer_ps = new LayerPoints[data->nz];
    for (INDEX k = 0; k < data->nz; k++) {
        ind_layer_ps[k].points[0].resize(data->nx * data->ny * 4);
        ind_layer_ps[k].points[1].resize(data->nx * data->ny * 4);
        ind_layer_ps[k].size[0] = 0;
        ind_layer_ps[k].size[1] = 0;
        for(INDEX j = 0; j < data->ny; j++) {
            for(INDEX i = 0; i < data->nx; i++) {
                auto &corner = data->GetCorner(k, j, i);
                for(INDEX layer = 0; layer < 2; layer++) { // z layer = 0 / 1
                    int is_added[2][2] = {{0, 0}, {0, 0}}; // 点是否已经添加
                    if(j > 0 && i > 0) {
                        auto &prev_corner = data->GetCorner(k, j-1, i-1);
                        if(isEqual(corner.zval[layer][0][0], prev_corner.zval[layer][1][1])) {
                            is_added[0][0] = 1;
                        }
                    }
                    if(j > 0) {
                        auto &prev_corner = data->GetCorner(k, j-1, i);
                        if(isEqual(corner.zval[layer][0][0], prev_corner.zval[layer][1][0])) {
                            is_added[0][0] = 1;
                        }
                        if(isEqual(corner.zval[layer][0][1], prev_corner.zval[layer][1][1])) {
                            is_added[0][1] = 1;
                        }
                    }
                    if(i > 0) {
                        auto &prev_corner = data->GetCorner(k, j, i-1);
                        if(isEqual(corner.zval[layer][0][0], prev_corner.zval[layer][0][1])) {
                            is_added[0][0] = 1;
                        }
                        if(isEqual(corner.zval[layer][1][0], prev_corner.zval[layer][1][1])) {
                            is_added[1][0] = 1;
                        }
                    }
                    for(int l = 0; l < 2; l++) {
                        for(int m = 0; m < 2; m++) {
                            if(is_added[l][m]) {
                                continue;
                            }
                            auto &p = corner.coord[layer][l][m];
                            ind_layer_ps[k].points[layer][ind_layer_ps[k].size[layer]] = p;
                            ind_layer_ps[k].size[layer]++;
                        }
                    }
                }
            }
        }
        ind_layer_ps[k].points[0].resize(ind_layer_ps[k].size[0]);
        ind_layer_ps[k].points[1].resize(ind_layer_ps[k].size[1]);
    }

    printf("finish adding independent layer points\n"); fflush(stdout);

    // Create a mesh from the points
    // mesh.vertices.assign_points(points, 3, data->nx * data->ny * data->nz); 

    layer_surfaces = new LayerSurface[data->nz - 1];
    for(int layer = 1; layer < data->nz; layer++) {
        layer_surfaces[layer - 1].pts.resize(data->nx * data->ny * 4);
        for(INDEX j = 0; j < data->ny; j++) {
            for(INDEX i = 0; i < data->nx; i++) {
                auto &corner = data->GetCorner(layer, j, i);
                layer_surfaces[layer - 1].pts[(j * data->nx + i) * 4 + 0] = corner.coord[0][0][0];
                layer_surfaces[layer - 1].pts[(j * data->nx + i) * 4 + 1] = corner.coord[0][0][1];
                layer_surfaces[layer - 1].pts[(j * data->nx + i) * 4 + 2] = corner.coord[0][1][0];
                layer_surfaces[layer - 1].pts[(j * data->nx + i) * 4 + 3] = corner.coord[0][1][1];
            }
        }
        for(int j = 0; j < data->ny; j++) {
            for(int i = 0; i < data->nx; i++) {
                auto &corner = data->GetCorner(layer, j, i);
                LayerSurface::Triangle t1 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 1, (j * data->nx + i)*4 + 3);
                LayerSurface::Triangle t2 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 2, (j * data->nx + i)*4 + 3);
                layer_surfaces[layer - 1].tris.push_back(t1);
                layer_surfaces[layer - 1].tris.push_back(t2);
                if(j > 0) {
                    auto &prev_corner = data->GetCorner(layer, j - 1, i);
                    if(!isEqual(corner.zval[0][0][0], prev_corner.zval[0][1][0]) && !isEqual(corner.zval[0][0][1], prev_corner.zval[0][1][1])) {
                        /// (j,i,0,0), (j,i,0,1), (j-1,i,1,0), (j-1,i,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 1, ((j-1) * data->nx + i)*4 + 3);
                        LayerSurface::Triangle t4 = LayerSurface::Triangle((j * data->nx + i)*4, ((j-1) * data->nx + i)*4 + 2, ((j-1) * data->nx + i)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                        layer_surfaces[layer - 1].tris.push_back(t4);
                    }
                    else if(!isEqual(corner.zval[0][0][1], prev_corner.zval[0][1][1])) {
                        /// (j,i,0,0), (j,i,0,1), (j-1,i,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 1, ((j-1) * data->nx + i)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                    }
                    else if(!isEqual(corner.zval[0][0][0], prev_corner.zval[0][1][0])) {
                        /// (j,i,0,0), (j-1,i,1,0), (j-1,i,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, ((j-1) * data->nx + i)*4 + 2, ((j-1) * data->nx + i)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                    }
                }
                if(i > 0) {
                    auto &prev_corner = data->GetCorner(layer, j, i - 1);
                    if(!isEqual(corner.zval[0][0][0], prev_corner.zval[0][0][1]) && !isEqual(corner.zval[0][1][0], prev_corner.zval[0][1][1])) {
                        /// (j,i,0,0), (j,i,1,0), (j,i-1,0,1), (j,i-1,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 2, (j * data->nx + i - 1)*4 + 3);
                        LayerSurface::Triangle t4 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i - 1)*4 + 1, (j * data->nx + i - 1)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                        layer_surfaces[layer - 1].tris.push_back(t4);
                    }
                    else if(!isEqual(corner.zval[0][1][0], prev_corner.zval[0][1][1])) {
                        /// (j,i,0,0), (j,i,1,0), (j,i-1,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i)*4 + 2, (j * data->nx + i - 1)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                    }
                    else if(!isEqual(corner.zval[0][0][0], prev_corner.zval[0][0][1])) {
                        /// (j,i,0,0), (j,i-1,0,1), (j,i-1,1,1)
                        LayerSurface::Triangle t3 = LayerSurface::Triangle((j * data->nx + i)*4, (j * data->nx + i - 1)*4 + 1, (j * data->nx + i - 1)*4 + 3);
                        layer_surfaces[layer - 1].tris.push_back(t3);
                    }
                }
            }
        }
        layer_surfaces[layer-1].UpdateTriRange();
        layer_surfaces[layer-1].BuildKDTree();
        // layer_surfaces[layer-1].SaveToVTK("layer_surfaces.vtk");
        printf("finish build layer surface %d\n", layer); fflush(stdout);
    }

    /// 3d Delaunay
    delaunay.periodic = false;
    // delaunay.periodic = true;
    delaunay.d = new GEO::PeriodicDelaunay3d(delaunay.periodic, 1.0);
    // delaunay.d = new GEO::ParallelDelaunay3d(3);
    if(!delaunay.periodic)
        delaunay.d->set_keeps_infinite(true);

    /// 2d Delaunay
    delaunay2d.d = new GEO::Delaunay2d();
}

double GeoGrid::GetPillarZvalAtLayer(int layer, Line pillar) {
    std::vector<LayerSurface::Triangle*> check_tris = layer_surfaces[layer - 1].kdtree.IntersectLine(pillar);
    for(int i = 0; i < check_tris.size(); i++) {
        Point3d p_res;
        int inct = geo3d::lineTriangleIntersection(pillar, check_tris[i]->p[0], check_tris[i]->p[1], check_tris[i]->p[2], p_res);
        if(inct > 0) {
            return p_res.z;
        }
    }
    printf("error : pillar not intersect with layer surface\n");
}

std::vector<Point3d> GeoGrid::GetMiddleGridPoints(Point3d p1, Point3d p2) {
    std::vector<Point3d> res;
    double dx = p2.x - p1.x;
    double dy = p2.y - p1.y;
    double dz = p2.z - p1.z;
    double dist = sqrt(dx * dx + dy * dy + dz * dz);
    int n = dist / grid_size;
    for(int i = 1; i < n; i++) {
        res.push_back(Point3d(p1.x + dx / dist * i * grid_size, 
            p1.y + dy / dist * i * grid_size,
            p1.z + dz / dist * i * grid_size));
    }
    return res;
}

void GeoGrid::GetSurfaceBorderFault(std::vector<Point3d> &pts, std::vector<INDEX> &bdrs) {
    // bottom surface is dense
    pts.reserve((data->ny + 1) * (data->nx + 1));
    fault_color = std::vector<int>((data->ny + 1) * (data->nx + 1), 0);
    // 颜色合并并查集结构
    std::vector<int> ufs_father((data->ny + 1) * (data->nx + 1), -1);
    for (int i = 0; i < (data->ny + 1) * (data->nx + 1); i++) {
        ufs_father[i] = i;
    }
    auto get_ufs_father = [&](int x) {
        while(x != ufs_father[x]) {
            x = ufs_father[x];
        }
        return x;
    };
    auto ufs_merge = [&](int x, int y) {
        int fx = get_ufs_father(x);
        int fy = get_ufs_father(y);
        if(fx != fy) {
            ufs_father[fx] = fy;
        }
    };

    for (INDEX j = 0; j < data->ny; j++) {
        for (INDEX i = 0; i < data->nx; i++) {
            pts.push_back(data->GetCorner(0, j, i).coord[0][0][0]);
        }
        pts.push_back(data->GetCorner(0, j, data->nx - 1).coord[0][0][1]);
    }
    for (INDEX i = 0; i < data->nx; i++) {
        pts.push_back(data->GetCorner(0, data->ny - 1, i).coord[0][1][0]);
    }
    pts.push_back(data->GetCorner(0, data->ny - 1, data->nx - 1).coord[0][1][1]);
    printf("points size: %d\n", pts.size());

    // add boundary lines, including four sides and fault lines

    // add side lines
    double x_side_length = data->range.rmax.x - data->range.rmin.x;
    double y_side_length = data->range.rmax.y - data->range.rmin.y;
    int x_side_num = x_side_length / grid_size * 2.0;
    int y_side_num = y_side_length / grid_size * 2.0;
    double x_side_size = x_side_length / x_side_num;
    double y_side_size = y_side_length / y_side_num;
    double z_val = data->GetCorner(0, 0, 0).zval[0][0][0];
    
    pts.push_back(Point3d(data->range.rmin.x, data->range.rmin.y, z_val));
    int start_idx = pts.size() - 1;
    for(int i = 0; i < x_side_num; i++) {
        double x_coord = (i == x_side_num - 1) ? data->range.rmax.x : data->range.rmin.x + (i + 1) * x_side_size;
        pts.push_back(Point3d(x_coord, data->range.rmin.y, z_val));
        bdrs.push_back(pts.size() - 1);
        bdrs.push_back(pts.size() - 2);
    }
    for(int i = 0; i < y_side_num; i++) {
        double y_coord = (i == y_side_num - 1) ? data->range.rmax.y : data->range.rmin.y + (i + 1) * y_side_size;
        pts.push_back(Point3d(data->range.rmax.x, y_coord, z_val));
        bdrs.push_back(pts.size() - 1);
        bdrs.push_back(pts.size() - 2);
    }
    for(int i = 0; i < x_side_num; i++) {
        double x_coord = (i == x_side_num - 1) ? 0 : data->range.rmax.x - (i + 1) * x_side_size;
        pts.push_back(Point3d(x_coord, data->range.rmax.y, z_val));
        bdrs.push_back(pts.size() - 1);
        bdrs.push_back(pts.size() - 2);
    }
    for(int i = 0; i < y_side_num - 1; i++) {
        double y_coord = (i == y_side_num - 1) ? 0 : data->range.rmax.y - (i + 1) * y_side_size;
        pts.push_back(Point3d(data->range.rmin.x, y_coord, z_val));
        bdrs.push_back(pts.size() - 1);
        bdrs.push_back(pts.size() - 2);
    }
    bdrs.push_back(pts.size() - 1);
    bdrs.push_back(start_idx);

    printf("side boundary line count: %d\n", bdrs.size() / 2);
    // for (int i = 0; i < bdrs.size(); i += 2) {
        // printf("boundary line at (%d %d)\n", bdrs[i], bdrs[i + 1]);
    // }

    // get fault points
    // int fault_pts = 0;
    // skip boundary points, only check z = 1 points
    // for (int k = 0; k < nz; k++) { // only consider the bottom surface
        for (int j = 0; j < data->ny + 1; j++) {
            for (int i = 0; i < data->nx + 1; i++) {
                if (j > 0 && j < data->ny && i > 0 && i < data->nx) { // common inner points
                    // if (!isEqual(data->GetCorner(0, j - 1, i - 1).zval[1][1][1], data->GetCorner(0, j, i - 1).zval[1][0][1]) ||
                    //     !isEqual(data->GetCorner(0, j - 1, i - 1).zval[1][1][1], data->GetCorner(0, j - 1, i).zval[1][1][0]) ||
                    //     !isEqual(data->GetCorner(0, j, i).zval[1][0][0], data->GetCorner(0, j - 1, i).zval[1][1][0]) ||
                    //     !isEqual(data->GetCorner(0, j, i).zval[1][0][0], data->GetCorner(0, j, i - 1).zval[1][0][1])) {
                    //     // printf("fault point at (%d, %d)\n", j, i);
                    //     data->isFault[j * (data->nx + 1) + i] = true;
                    //     fault_pts++;
                    //     continue;
                    // }
                    if (!isEqual(data->GetCorner(0, j - 1, i - 1).zval[1][1][1], data->GetCorner(0, j, i - 1).zval[1][0][1])) {
                        data->isFault[j * (data->nx + 1) + i] = true;
                        data->isFault[j * (data->nx + 1) + i - 1] = true;
                    }
                    if (!isEqual(data->GetCorner(0, j - 1, i - 1).zval[1][1][1], data->GetCorner(0, j - 1, i).zval[1][1][0])) {
                        data->isFault[j * (data->nx + 1) + i] = true;
                        data->isFault[(j - 1) * (data->nx + 1) + i] = true;
                    }
                    if (!isEqual(data->GetCorner(0, j, i).zval[1][0][0], data->GetCorner(0, j - 1, i).zval[1][1][0])) {
                        data->isFault[j * (data->nx + 1) + i] = true;
                        data->isFault[j * (data->nx + 1) + i + 1] = true;
                    }
                    if (!isEqual(data->GetCorner(0, j, i).zval[1][0][0], data->GetCorner(0, j, i - 1).zval[1][0][1])) {
                        data->isFault[j * (data->nx + 1) + i] = true;
                        data->isFault[(j + 1) * (data->nx + 1) + i] = true;
                    }
                }
                else {
                    // only have inner fault, skip
                }
            }
        }
    // }
    // printf("fault points count: %d\n", fault_pts);
    // Warning: the st/ed fault points are not correct

    // 断层点连边
    std::vector<std::vector<int>> fault_edges;
    // 是否为关键断层点
    std::vector<bool> is_curcial((data->ny + 1) * (data->nx + 1), false);
    fault_edges.resize((data->ny + 1) * (data->nx + 1));
    // 等待处理的断层链代表的边
    std::queue<std::pair<int, int> > path_que;

    // add fault lines
    int fault_lines = 0;
    for (int j = 0; j < data->ny + 1; j++) {
        for (int i = 0; i < data->nx + 1; i++) {
            if (data->isFault[j * (data->nx + 1) + i]) {
                if (data->isFault[j * (data->nx + 1) + i + 1]) {
                    fault_edges[j * (data->nx + 1) + i].push_back(j * (data->nx + 1) + i + 1);
                    fault_edges[j * (data->nx + 1) + i + 1].push_back(j * (data->nx + 1) + i);

                    // bdrs.push_back(j * (data->nx + 1) + i);
                    // bdrs.push_back(j * (data->nx + 1) + i + 1);

                    fault_lines++;
                }
                if (data->isFault[(j + 1) * (data->nx + 1) + i]) {
                    // printf("fault line at (%d, %d) and (%d, %d)\n", j, i, j + 1, i);
                    fault_edges[j * (data->nx + 1) + i].push_back((j + 1) * (data->nx + 1) + i);
                    fault_edges[(j + 1) * (data->nx + 1) + i].push_back(j * (data->nx + 1) + i);

                    // bdrs.push_back(j * (data->nx + 1) + i);
                    // bdrs.push_back((j + 1) * (data->nx + 1) + i);

                    fault_lines++;
                }
            }
        }
    }
    for (int j = 0; j < data->ny + 1; j++) {
        for (int i = 0; i < data->nx + 1; i++) {
            if (data->isFault[j * (data->nx + 1) + i]) {
                // with degree 1 or degree > 2
                if (fault_edges[j * (data->nx + 1) + i].size() == 1 || fault_edges[j * (data->nx + 1) + i].size() > 2) {
                    is_curcial[j * (data->nx + 1) + i] = true;
                    if(fault_edges[j * (data->nx + 1) + i].size() == 1 && path_que.empty()) {
                        path_que.push(std::make_pair(j * (data->nx + 1) + i, fault_edges[j * (data->nx + 1) + i][0]));
                    }
                }
                else {
                    // degree = 2
                    double cos_angle = geo3d::cosAngle(pts[j * (data->nx + 1) + i], pts[fault_edges[j * (data->nx + 1) + i][0]], pts[fault_edges[j * (data->nx + 1) + i][1]]);
                    if(cos_angle >= cos_fault_angle_threshold) {
                        is_curcial[j * (data->nx + 1) + i] = true;
                        ufs_merge(j * (data->nx + 1) + i, fault_edges[j * (data->nx + 1) + i][0]);
                        ufs_merge(j * (data->nx + 1) + i, fault_edges[j * (data->nx + 1) + i][1]);
                    }
                }
            }
        }
    }

    for(int i = 0; i < (data->ny + 1) * (data->nx + 1); i++) {
        if (fault_edges[i].size() > 2) {
            for(int j = 0; j < fault_edges[i].size(); j++) {
                for(int k = j + 1; k < fault_edges[i].size(); k++) {
                    double cos_angle = geo3d::cosAngle(pts[i], pts[fault_edges[i][j]], pts[fault_edges[i][k]]);
                    // 夹角较平，属于同一断层
                    if (cos_angle < cos_fault_angle_threshold) {
                        ufs_merge(i, fault_edges[i][j]);
                        ufs_merge(i, fault_edges[i][k]);
                    }
                }
            }
        }
        else if(fault_edges[i].size() == 1) {
            ufs_merge(i, fault_edges[i][0]);
        }
    }

    while(!path_que.empty()) {
        auto path = path_que.front();
        path_que.pop();
        int st = path.first;
        int ed = path.second;
        // get the path of st, ed
        int prev = st;
        int curr = ed;
        std::vector<int> path_list;
        path_list.push_back(st);
        path_list.push_back(ed);
        while(!is_curcial[curr]) {
            int next = -1;
            for(int i = 0; i < fault_edges[curr].size(); i++) {
                if(fault_edges[curr][i] != prev) {
                    next = fault_edges[curr][i];
                    break;
                }
            }
            prev = curr;
            curr = next;
            path_list.push_back(curr);
        }
        // add path of next curcial point to queue
        for(int i = 0; i < fault_edges[curr].size(); i++) {
            if(fault_edges[curr][i] != prev) {
                path_que.push(std::make_pair(curr, fault_edges[curr][i]));
            }
        }
        // 并查集合并路径中间节点
        for(int i = 2; i < path_list.size() - 1; i++) {
            ufs_merge(path_list[i], path_list[i - 1]);
        }
        // remove redundant points from queue
        int remove_index;
        do {
            remove_index = -1;
            double new_dist = 0;
            for(int i = 1; i < path_list.size() - 1; i++) {
                int edge_len0 = geo3d::distance(pts[path_list[i - 1]], pts[path_list[i]]);
                int edge_len1 = geo3d::distance(pts[path_list[i]], pts[path_list[i + 1]]);
                // one of the edge is too short
                if(edge_len0 < grid_size * 0.5 || edge_len1 < grid_size * 0.5) {
                    double cos_angle = geo3d::cosAngle(pts[path_list[i]], pts[path_list[i - 1]], pts[path_list[i + 1]]);
                    if(cos_angle >= cos_fault_angle_threshold) {
                        // 夹角较大，不可删除
                        continue;
                    }
                    // 考虑删除
                    if(remove_index != -1 && new_dist < edge_len0 + edge_len1) {
                        // 删除上一个点较好
                        continue;
                    }
                    remove_index = i;
                    new_dist = edge_len0 + edge_len1;
                }
            }
            if(remove_index != -1) {
                path_list.erase(path_list.begin() + remove_index);
            }
        }while(remove_index != -1);
        // add path to boundary
        for(int i = 0; i < path_list.size() - 1; i++) {
            bdrs.push_back(path_list[i]);
            bdrs.push_back(path_list[i + 1]);
        }
    }

    // 不同断层数量
    int fault_num = 0;
    for(int i = 0; i < (data->ny + 1) * (data->nx + 1); i++) {
        if(ufs_father[i] == i && data->isFault[i]) {
            fault_num++;
            for(int j = 0; j < (data->ny + 1) * (data->nx + 1); j++) {
                if(get_ufs_father(j) == i) {
                    fault_color[j] = fault_num;
                }
            }
        }
    }
    printf("number of different faults : %d\n", fault_num);

    printf("boundary line count: %d\n", bdrs.size() / 2);    

    std::ofstream vtkFile("surface.vtk");

    // Write the VTK file header
    vtkFile << "# vtk DataFile Version 2.0\n";
    vtkFile << "Grid Legacy VTK File\n";
    vtkFile << "ASCII\n";
    vtkFile << "DATASET POLYDATA\n";
    vtkFile << "\n";

    // Write the grid points to the VTK file
    vtkFile << "POINTS " << pts.size() << " float\n";
    for(int i = 0; i < pts.size(); i++) {
        vtkFile << pts[i].x << " " << pts[i].y << " " << pts[i].z << "\n";
    }
    vtkFile << "\n";
    // Write the grid cells to the VTK file
    vtkFile << "LINES " << bdrs.size() / 2 << " " << bdrs.size() / 2 * 3 << "\n";
    for(int i = 0; i < bdrs.size() / 2; i++) {
        vtkFile << 2 << " " << bdrs[2*i] << " " << bdrs[2*i + 1] << "\n";
    }
    vtkFile << "\n";
    vtkFile << "POINT_DATA " << pts.size() << "\n";
    vtkFile << "SCALARS isInterface int 1\nLOOKUP_TABLE default\n";
    for(int i = 0; i < (data->ny + 1) * (data->nx + 1); i++) {
        int val = 0;
        if(data->isFault[i]) val = 1;
        if(is_curcial[i]) val = 2;
        vtkFile << val << "\n";
    }
    for(int i = (data->ny + 1) * (data->nx + 1); i < pts.size(); i++) {
        vtkFile << 0 << "\n";
    }
    // Close the VTK file
    vtkFile.close();
}

void GeoGrid::SetSurfaceByBottomSurface(const std::vector<Point3d> &pts, const std::vector<INDEX> &bdrs) {
    /// 构建 2d Delaunay
    std::vector<INDEX> indices = bdrs;
    std::sort(indices.begin(), indices.end());
    auto it = std::unique(indices.begin(), indices.end());
    indices.erase(it, indices.end());
    delaunay2d.n_points = indices.size();
    delaunay2d.flat_points.resize(2 * delaunay2d.n_points);
    for(int i = 0; i < delaunay2d.n_points; i++) {
        delaunay2d.flat_points[2*i] = pts[indices[i]].x;
        delaunay2d.flat_points[2*i+1] = pts[indices[i]].y;
    }

    printf("delaunay2d set_vertices\n"); fflush(stdout);
    delaunay2d.d->set_vertices(delaunay2d.n_points, delaunay2d.flat_points.data());

    // printf("delaunay2d init SaveToVTK\n"); fflush(stdout);
    // delaunay2d.SaveToVTK("../../output/delaunay2d_0.vtk");
    // printf("delaunay2d finish SaveToVTK\n"); fflush(stdout);

    int iter = 0;
    // 二维 Delaunay 迭代网格细化
    while(iter < 42) {
        iter++;
        printf("refinement iter %d\n", iter);
        int inc_point = 0;
        std::vector<double> new_points;
        inc_point = 0;
        int n_trangles = delaunay2d.d->nb_cells();
        for(int i = 0; i < n_trangles; i++) {
            const GEO::signed_index_t* cell = delaunay2d.d->cell_to_v() + 3*i;
            geo2d::Point2d p[3] = {delaunay2d.getPointByIndex(cell[0]), delaunay2d.getPointByIndex(cell[1]), delaunay2d.getPointByIndex(cell[2])};
            double radius;
            geo2d::Point2d center;
            geo2d::getTriangleCircumcenter(p[0], p[1], p[2], center, radius);
            if(radius >= grid_size) {
                bool flag = false;
                for(int j = 0; j < inc_point; j++) {
                    double dist = sqrt((new_points[j * 2] - center.x) * (new_points[j * 2] - center.x) + (new_points[j * 2 + 1] - center.y) * (new_points[j * 2 + 1] - center.y));
                    if(dist < grid_size) {
                        flag = true;
                        break;
                    }
                }
                for(int j = 0; j < bdrs.size() / 2; j++) {
                    geo2d::Point2d p1 = geo2d::Point2d(pts[bdrs[j * 2]].x, pts[bdrs[j * 2]].y);
                    geo2d::Point2d p2 = geo2d::Point2d(pts[bdrs[j * 2 + 1]].x, pts[bdrs[j * 2 + 1]].y);
                    if(geo2d::InCircle(p1, p2, center)) {
                        flag = true;
                        break;
                    }
                }
                if(!flag) {
                    new_points.push_back(center.x);
                    new_points.push_back(center.y);
                    inc_point++;
                }
            }
        }
        if(inc_point == 0) {
            printf("refinement finish at iter %d\n", iter);
            break;
        }
        printf("refinement iter %d add %d points\n", iter, inc_point);
        delaunay2d.n_points += inc_point;
        delaunay2d.flat_points.insert(delaunay2d.flat_points.end(), new_points.begin(), new_points.end());
        delaunay2d.d->set_vertices(delaunay2d.n_points, delaunay2d.flat_points.data());

        // std::string filename = "../../output/delaunay2d_" + std::to_string(iter) + ".vtk";
        // delaunay2d.SaveToVTK(filename.c_str());
    }

    /// 构建 3d 点集合
    double pts_zval = pts[0].z;
    double top_level_zval = data->GetCorner(0, 0, 0).zval[0][0][0];
    double bottom_level_zval = data->GetCorner(data->nz - 1, 0, 0).zval[1][1][1];

    std::vector<Point3d>* pillar_points = new std::vector<Point3d>[delaunay2d.n_points];
    std::vector<int>* pillar_point_indices = new std::vector<int>[delaunay2d.n_points];  
    Line* pillar_lines = new Line[delaunay2d.n_points];
    int index_count = 0;

    for(int i = 0; i < delaunay2d.n_points; i++) {
        std::vector<double> zvals;
        if(i < indices.size() && data->isFault[indices[i]]) {
            int idx = indices[i];
            // fault pillar, add more than one points
            int y_idx = idx / (data->nx + 1);
            int x_idx = idx % (data->nx + 1);
            if(y_idx <= 0 || y_idx >= data->ny || x_idx <= 0 || x_idx >= data->nx) {
                printf("invalid fault point on border\n"); exit(-1);
            }
            INDEX px, py;
            data->InterpolateNewPillar(delaunay2d.flat_points[2*i], delaunay2d.flat_points[2*i+1], pts_zval, pillar_lines[i], px, py);
            zvals.push_back(top_level_zval);
            zvals.push_back(bottom_level_zval);
            for(int layer = 1; layer < data->nz; layer++) {
                zvals.push_back(data->GetCorner(layer, y_idx, x_idx).zval[0][0][0]);
                zvals.push_back(data->GetCorner(layer, y_idx, x_idx-1).zval[0][0][1]);
                zvals.push_back(data->GetCorner(layer, y_idx-1, x_idx).zval[0][1][0]);
                zvals.push_back(data->GetCorner(layer, y_idx-1, x_idx-1).zval[0][1][1]);
            }
        }
        else {
            INDEX px, py;
            data->InterpolateNewPillar(delaunay2d.flat_points[2*i], delaunay2d.flat_points[2*i+1], pts_zval, pillar_lines[i], px, py);
            zvals.push_back(top_level_zval);
            zvals.push_back(bottom_level_zval);
            for(int layer = 1; layer < data->nz; layer++) {
                zvals.push_back(GetPillarZvalAtLayer(layer, pillar_lines[i]));
            }
        }
        std::sort(zvals.begin(), zvals.end());
        std::vector<Point3d> pillar_pts;
        for(int j = 0; j < zvals.size(); j++) {
            if(j == 0 || !isEqual(zvals[j], zvals[j-1])) {
                geo2d::Point2d p_curr = data->GetPillarAtZval(pillar_lines[i], zvals[j]);
                pillar_pts.push_back(Point3d(p_curr.x, p_curr.y, zvals[j]));
                pillar_points[i].push_back(Point3d(p_curr.x, p_curr.y, zvals[j]));
                pillar_point_indices[i].push_back(index_count++);
                if(i < indices.size() && data->isFault[indices[i]]) {
                    // 断层点
                    // mesh_points_type.push_back(PointType::PT_FAULT);
                    mesh_points_type.push_back(fault_color[indices[i]]);
                }
                else if(zvals[j] != top_level_zval && zvals[j] != bottom_level_zval) {
                    // 地层点
                    // mesh_points_type.push_back(PointType::PT_LAYER);
                    mesh_points_type.push_back(10 + j + 1);
                }
                // }
                else {
                    // 其他点
                    // mesh_points_type.push_back(PointType::PT_OTHER);
                    mesh_points_type.push_back(0);
                }
            }
        }
        for(int j = 1; j < pillar_pts.size(); j++) {
            std::vector<Point3d> middle_points = GetMiddleGridPoints(pillar_pts[j - 1], pillar_pts[j]);
            for(int k = 0; k < middle_points.size(); k++) {
                pillar_points[i].push_back(middle_points[k]);
                pillar_point_indices[i].push_back(index_count++);
                if(i < indices.size() && data->isFault[indices[i]]) {
                    // 断层点
                    // mesh_points_type.push_back(PointType::PT_FAULT);
                    mesh_points_type.push_back(fault_color[indices[i]]);
                }
                else {
                    // 其他点
                    // mesh_points_type.push_back(PointType::PT_OTHER);
                    mesh_points_type.push_back(0);
                }
            }
        }        
    }
    for(int i = 0; i < delaunay2d.n_points; i++) {
        for(int j = 0; j < pillar_points[i].size(); j++) {
            mesh_points.push_back(pillar_points[i][j]);
        }
    }
    printf("add a total of %d points\n", index_count);

    if(strategy == 0) {
        // 使用所有点直接划分
        std::vector<Point3d> surface_pts;
        for(int i = 0; i < delaunay2d.n_points; i++) {
            for(int j = 0; j < pillar_points[i].size(); j++) {
                surface_pts.push_back(pillar_points[i][j]);
            }
        }

        printf("begin delaunay 3d calculation\n"); fflush(stdout);

        delaunay.flat_points.resize(3 * surface_pts.size());
        for(int i = 0; i < surface_pts.size(); i++) {
            delaunay.flat_points[3*i] = surface_pts[i].x;
            delaunay.flat_points[3*i+1] = surface_pts[i].y;
            delaunay.flat_points[3*i+2] = surface_pts[i].z;
        }
        delaunay.d->set_vertices(surface_pts.size(), delaunay.flat_points.data());
        delaunay.d->compute();

        printf("after delaunay 3d calculation\n"); fflush(stdout);

        int n_cells = delaunay.d->nb_cells();
        for(int i = 0; i < n_cells; i++) {
            GEO::signed_index_t* cell = (GEO::signed_index_t*)(delaunay.d->cell_to_v()) + 4*i;
            bool flag = 0;
            for(int j = 0; j < 4; j++) {
                int idx = cell[j];
                if(idx < 0) {
                    flag = 1;
                }
            }
            if(flag) continue;
            for(int j = 0; j < 4; j++) {
                tetrahedron_cells.push_back(cell[j]);
            }
        }
        printf("point number : %d\n", surface_pts.size());
        printf("tetrahedron cell number : %d\n", tetrahedron_cells.size() / 4);
    }
    else if(strategy == 1) {
        // 对于每层分别划分
        std::vector<Point3d> surface_pts;
        std::vector<int> surface_pt_indices;
        std::vector<PointType> surface_pts_type;
        for(int layer = 0; layer < data->nz; layer++) {
            // perform delaunay 3d for layer
            printf("perform delaunay 3d for layer %d\n", layer);
            surface_pts.clear();
            surface_pts_type.clear();
            surface_pt_indices.clear();
            for(int i = 0; i < delaunay2d.n_points; i++) {
                std::vector<double> zvals;
                if(i < indices.size() && data->isFault[indices[i]]) {
                    int idx = indices[i];
                    // fault, add more than one points
                    int y_idx = idx / (data->nx + 1);
                    int x_idx = idx % (data->nx + 1);
                    if(y_idx <= 0 || y_idx >= data->ny || x_idx <= 0 || x_idx >= data->nx) {
                        printf("invalid fault point on border\n"); exit(-1);
                    }
                    if(layer == 0) {
                        zvals.push_back(data->GetCorner(0, y_idx, x_idx).zval[1][0][0]);
                        zvals.push_back(data->GetCorner(0, y_idx, x_idx-1).zval[1][0][1]);
                        zvals.push_back(data->GetCorner(0, y_idx-1, x_idx).zval[1][1][0]);
                        zvals.push_back(data->GetCorner(0, y_idx-1, x_idx-1).zval[1][1][1]);
                        zvals.push_back(top_level_zval);
                    }
                    else if(layer == data->nz - 1) {
                        zvals.push_back(data->GetCorner(data->nz-1, y_idx, x_idx).zval[0][0][0]);
                        zvals.push_back(data->GetCorner(data->nz-1, y_idx, x_idx-1).zval[0][0][1]);
                        zvals.push_back(data->GetCorner(data->nz-1, y_idx-1, x_idx).zval[0][1][0]);
                        zvals.push_back(data->GetCorner(data->nz-1, y_idx-1, x_idx-1).zval[0][1][1]);
                        zvals.push_back(bottom_level_zval);
                    }
                    else {
                        zvals.push_back(data->GetCorner(layer-1, y_idx, x_idx).zval[1][0][0]);
                        zvals.push_back(data->GetCorner(layer-1, y_idx, x_idx-1).zval[1][0][1]);
                        zvals.push_back(data->GetCorner(layer-1, y_idx-1, x_idx).zval[1][1][0]);
                        zvals.push_back(data->GetCorner(layer-1, y_idx-1, x_idx-1).zval[1][1][1]);
                        zvals.push_back(data->GetCorner(layer, y_idx, x_idx).zval[1][0][0]);
                        zvals.push_back(data->GetCorner(layer, y_idx, x_idx-1).zval[1][0][1]);
                        zvals.push_back(data->GetCorner(layer, y_idx-1, x_idx).zval[1][1][0]);
                        zvals.push_back(data->GetCorner(layer, y_idx-1, x_idx-1).zval[1][1][1]);
                    }
                    std::sort(zvals.begin(), zvals.end());
                    for(int j = 0; j < pillar_points[i].size(); j++) {
                        if(pillar_points[i][j].z >= zvals[0] && pillar_points[i][j].z <= zvals[zvals.size() - 1]) {
                            surface_pts.push_back(pillar_points[i][j]);
                            surface_pt_indices.push_back(pillar_point_indices[i][j]);
                            if(pillar_points[i][j].z == top_level_zval) {
                                surface_pts_type.push_back(PointType::PT_UPPER_SURFACE);
                            }
                            else if(pillar_points[i][j].z == bottom_level_zval) {
                                surface_pts_type.push_back(PointType::PT_LOWER_SURFACE);
                            }
                            else {
                                surface_pts_type.push_back(PointType::PT_FAULT);
                            }
                        }
                    }
                }
                else {
                    double zval_up;
                    double zval_down;
                    if(layer == 0) {
                        zval_up = top_level_zval;
                        zval_down = GetPillarZvalAtLayer(1, pillar_lines[i]);
                    }
                    else if(layer == data->nz - 1) {
                        zval_up = GetPillarZvalAtLayer(data->nz - 1, pillar_lines[i]);
                        zval_down = bottom_level_zval;
                    }
                    else {
                        zval_up = GetPillarZvalAtLayer(layer, pillar_lines[i]);
                        zval_down = GetPillarZvalAtLayer(layer + 1, pillar_lines[i]);
                    }
                    for(int j = 0; j < pillar_points[i].size(); j++) {
                        if((pillar_points[i][j].z >= zval_down && pillar_points[i][j].z <= zval_up) ||
                            (pillar_points[i][j].z >= zval_up && pillar_points[i][j].z <= zval_down)) {
                            surface_pts.push_back(pillar_points[i][j]);
                            surface_pt_indices.push_back(pillar_point_indices[i][j]);
                            if(pillar_points[i][j].z == zval_up) {
                                surface_pts_type.push_back(PointType::PT_UPPER_SURFACE);
                            }
                            else if(pillar_points[i][j].z == zval_down) {
                                surface_pts_type.push_back(PointType::PT_LOWER_SURFACE);
                            }
                            else {
                                surface_pts_type.push_back(PointType::PT_INNER);
                            }
                        }
                    }
                }
            }

            printf("begin delaunay 3d calculation\n"); fflush(stdout);

            delaunay.flat_points.resize(3 * surface_pts.size());
            for(int i = 0; i < surface_pts.size(); i++) {
                delaunay.flat_points[3*i] = surface_pts[i].x;
                delaunay.flat_points[3*i+1] = surface_pts[i].y;
                delaunay.flat_points[3*i+2] = surface_pts[i].z;
            }
            delaunay.d->set_vertices(surface_pts.size(), delaunay.flat_points.data());
            delaunay.d->compute();

            printf("after delaunay 3d compute\n"); fflush(stdout);

            int n_cells = delaunay.d->nb_cells();
            const double* vertices = delaunay.d->vertices_ptr();
            int set_invalid_cells = 0;
            for(int i = 0; i < n_cells; i++) {
                GEO::signed_index_t* cell = (GEO::signed_index_t*)(delaunay.d->cell_to_v()) + 4*i;
                std::vector<Point3d> ps;
                bool flag = 0;
                for(int j = 0; j < 4; j++) {
                    int idx = cell[j];
                    if(idx < 0) {
                        flag = 1;
                    } else {
                        ps.push_back(Point3d(vertices[idx * 3], vertices[idx * 3 + 1], vertices[idx * 3 + 2]));
                    }
                }
                if(flag) continue;
                // 去除所有点都在上表面或下表面的四面体
                if((surface_pts_type[cell[0]] == PT_UPPER_SURFACE || surface_pts_type[cell[0]] == PT_LOWER_SURFACE) &&
                    surface_pts_type[cell[1]] == surface_pts_type[cell[0]] && surface_pts_type[cell[2]] == surface_pts_type[cell[0]] && surface_pts_type[cell[3]] == surface_pts_type[cell[0]]) {
                    cell[0] = -1;
                    set_invalid_cells++;
                    continue;
                }

                /// 四面体中心
                Point3d p_incenter = geo3d::computeIncenter(ps[0], ps[1], ps[2], ps[3]);
                /// 四面体重心
                Point3d p_masscenter = geo3d::computeMasscenter(ps[0], ps[1], ps[2], ps[3]);

                if(layer != data->nz - 1) {
                    // 与本层上表面判断交点数
                    // 去除距离表面过近的四面体
                    p_incenter.z += grid_size * 0.01;
                    // 四面体中心向上的射线与地层的交点数
                    int insert_count = layer_surfaces[layer].checkInsertWithUpLine(p_incenter);
                    if(insert_count % 2 == 0) {
                        // 四面体中心在底层外部
                        cell[0] = -1;
                        set_invalid_cells++;
                        continue;
                    }
                    // 去除距离表面过近的四面体
                    p_masscenter.z += grid_size * 0.01;
                    // 四面体中心向上的射线与地层的交点数
                    insert_count = layer_surfaces[layer].checkInsertWithUpLine(p_masscenter);
                    if(insert_count % 2 == 0) {
                        // 四面体中心在底层外部
                        cell[0] = -1;
                        set_invalid_cells++;
                        continue;
                    }
                }
                if(layer != 0) {
                    // 与本层下表面判断交点数
                    // 去除距离表面过近的四面体
                    p_incenter.z -= grid_size * 0.01;
                    // 四面体中心向下的射线与地层的交点数
                    int insert_count = layer_surfaces[layer - 1].checkInsertWithDownLine(p_incenter);
                    if(insert_count % 2 == 0) {
                        // 四面体中心在底层外部
                        cell[0] = -1;
                        set_invalid_cells++;
                        continue;
                    }
                    // 去除距离表面过近的四面体
                    p_masscenter.z -= grid_size * 0.01;
                    // 四面体中心向下的射线与地层的交点数
                    insert_count = layer_surfaces[layer - 1].checkInsertWithDownLine(p_masscenter);
                    if(insert_count % 2 == 0) {
                        // 四面体中心在底层外部
                        cell[0] = -1;
                        set_invalid_cells++;
                        continue;
                    }
                }

                for(int j = 0; j < 4; j++) {
                    tetrahedron_cells.push_back(surface_pt_indices[cell[j]]);
                }
            }
            printf("set %d cells to invalid\n", set_invalid_cells);
            std::string output_file = data->output_path + "/delaunay_layer" + std::to_string(layer) + ".vtk"; 
            printf("saving to %s\n", output_file.c_str());
            delaunay.SaveToVTK(output_file);
        }
    }

    // printf("saving full mesh\n");
    // Mesh3D mesh;
    // for(int i = 0; i < delaunay2d.n_points; i++) {
    //     for(int j = 0; j < pillar_points[i].size(); j++) {
    //         mesh.points.push_back(pillar_points[i][j]);
    //     }
    // }
    // mesh.tetrahedron_cells = tetrahedron_cells;
    // mesh.young.resize(mesh.tetrahedron_cells.size() / 4, 0.0f);
    // mesh.poisson.resize(mesh.tetrahedron_cells.size() / 4, 0.0f);
    // mesh.rockden.resize(mesh.tetrahedron_cells.size() / 4, 0.0f);
    // mesh.poro.resize(mesh.tetrahedron_cells.size() / 4, 0.0f);
    // mesh.WriteVTK("../../output/delaunay.vtk");

    delete[] pillar_points;
    delete[] pillar_point_indices;
    delete[] pillar_lines;
}

void GeoGrid::Expand() {
    /// top 方向数值较小
    assert(expand_top < expand_bottom);
    is_expand = true;

    mesh_points_expand = mesh_points;
    tetrahedron_cells_expand = tetrahedron_cells;

    expand_info.inner_rmin = data->range.rmin;
    expand_info.inner_rmax = data->range.rmax;
    
    Point3d inner_size = expand_info.inner_rmax - expand_info.inner_rmin;
    expand_info.outer_rmin.z = expand_info.expand_top;
    expand_info.outer_rmax.z = expand_info.expand_bottom;
    expand_info.outer_rmin.x = expand_info.inner_rmin.x - inner_size.x * expand_info.expand_X_ratio;
    expand_info.outer_rmax.x = expand_info.inner_rmax.x + inner_size.x * expand_info.expand_X_ratio;
    expand_info.outer_rmin.y = expand_info.inner_rmin.y - inner_size.y * expand_info.expand_Y_ratio;
    expand_info.outer_rmax.y = expand_info.inner_rmax.y + inner_size.y * expand_info.expand_Y_ratio;

    int x_expand_layers = inner_size.x * expand_info.expand_X_ratio / grid_size;
    double x_expand_grid_size = inner_size.x * expand_info.expand_X_ratio / x_expand_layers;
    int y_expand_layers = inner_size.y * expand_info.expand_Y_ratio / grid_size;
    double y_expand_grid_size = inner_size.y * expand_info.expand_Y_ratio / y_expand_layers;
    int top_expand_layers = (expand_info.inner_rmin.z - expand_info.outer_rmin.z) / grid_size;
    double top_expand_grid_size = (expand_info.inner_rmin.z - expand_info.outer_rmin.z) / top_expand_layers;
    int bottom_expand_layers = (expand_info.outer_rmax.z - expand_info.inner_rmax.z) / grid_size;
    double bottom_expand_grid_size = (expand_info.outer_rmax.z - expand_info.inner_rmax.z) / bottom_expand_layers;

    int left_expand_points = 0;
    int left_points = 0;

    // 补充前后左右四块中的点
    for(int i = 0; i < mesh_points.size(); i++) {
        Point3d p = mesh_points[i];
        if(p.x == expand_info.inner_rmin.x) {
            // 左侧点
            left_points++;
            for(int j = 0; j < x_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x - x_expand_grid_size * (j + 1), p.y, p.z));
                left_expand_points++;
            }
        }
        if(p.x == expand_info.inner_rmax.x) {
            // 右侧点
            for(int j = 0; j < x_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x + x_expand_grid_size * (j + 1), p.y, p.z));
            }
        }
        if(p.y == expand_info.inner_rmin.y) {
            // 前侧点
            for(int j = 0; j < y_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x, p.y - y_expand_grid_size * (j + 1), p.z));
            }
        }
        if(p.y == expand_info.inner_rmax.y) {
            // 后侧点
            for(int j = 0; j < y_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x, p.y + y_expand_grid_size * (j + 1), p.z));
            }
        }
    }

    // 补充四个斜方向中的点
    for(int i = 0; i < mesh_points.size(); i++) {
        Point3d p = mesh_points[i];
        if(p.x == expand_info.inner_rmin.x && p.y == expand_info.inner_rmin.y) {
            // 左前角点
            for(int j = 0; j < x_expand_layers; j++) {
                for(int k = 0; k < y_expand_layers; k++) {
                    mesh_points_expand.push_back(Point3d(p.x - x_expand_grid_size * (j + 1), p.y - y_expand_grid_size * (k + 1), p.z));
                }
            }
        }
        if(p.x == expand_info.inner_rmin.x && p.y == expand_info.inner_rmax.y) {
            // 左后角点
            for(int j = 0; j < x_expand_layers; j++) {
                for(int k = 0; k < y_expand_layers; k++) {
                    mesh_points_expand.push_back(Point3d(p.x - x_expand_grid_size * (j + 1), p.y + y_expand_grid_size * (k + 1), p.z));
                }
            }
        }
        if(p.x == expand_info.inner_rmax.x && p.y == expand_info.inner_rmin.y) {
            // 右前角点
            for(int j = 0; j < x_expand_layers; j++) {
                for(int k = 0; k < y_expand_layers; k++) {
                    mesh_points_expand.push_back(Point3d(p.x + x_expand_grid_size * (j + 1), p.y - y_expand_grid_size * (k + 1), p.z));
                }
            }
        }
        if(p.x == expand_info.inner_rmax.x && p.y == expand_info.inner_rmax.y) {
            // 右后角点
            for(int j = 0; j < x_expand_layers; j++) {
                for(int k = 0; k < y_expand_layers; k++) {
                    mesh_points_expand.push_back(Point3d(p.x + x_expand_grid_size * (j + 1), p.y + y_expand_grid_size * (k + 1), p.z));
                }
            }
        }
    }
    
    int cur_tot_points = mesh_points_expand.size();
    // 补充上方下方点
    for(int i = 0; i < cur_tot_points; i++) {
        Point3d p = mesh_points_expand[i];
        if(p.z == expand_info.inner_rmin.z) {
            // 上方点
            for(int j = 0; j < top_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x, p.y, p.z - top_expand_grid_size * (j + 1)));
            }
        }
        if(fabs(p.z - expand_info.inner_rmax.z) < 1e-6) {
            // 下方点
            for(int j = 0; j < bottom_expand_layers; j++) {
                mesh_points_expand.push_back(Point3d(p.x, p.y, p.z + bottom_expand_grid_size * (j + 1)));
            }
        }
    }
    printf("expanded point number : %d\n", mesh_points_expand.size());

    delaunay.flat_points.resize(3 * mesh_points_expand.size());
    for(int i = 0; i < mesh_points_expand.size(); i++) {
        delaunay.flat_points[3*i] = mesh_points_expand[i].x;
        delaunay.flat_points[3*i+1] = mesh_points_expand[i].y;
        delaunay.flat_points[3*i+2] = mesh_points_expand[i].z;
    }

    printf("begin delaunay 3d expand calculation\n"); fflush(stdout);

    delaunay.d->set_vertices(mesh_points_expand.size(), delaunay.flat_points.data());
    delaunay.d->compute();
    
    printf("after delaunay 3d expand calculation\n"); fflush(stdout);

    // std::string output_file = data->output_path + "/delaunay_expand.vtk"; 
    // printf("saving to %s\n", output_file.c_str());
    // delaunay.SaveToVTK(output_file); 

    int n_cells = delaunay.d->nb_cells();
    for(int i = 0; i < n_cells; i++) {
        GEO::signed_index_t* cell = (GEO::signed_index_t*)(delaunay.d->cell_to_v()) + 4*i;
        bool flag = 0;
        for(int j = 0; j < 4; j++) {
            int idx = cell[j];
            if(idx < 0) {
                flag = 1;
            }
        }
        if(flag) continue;
        Point3d direction = expand_info.GetTetrahedronPosition(mesh_points_expand[cell[0]], mesh_points_expand[cell[1]], mesh_points_expand[cell[2]], mesh_points_expand[cell[3]]);
        if(direction.x == 0 && direction.y == 0 && direction.z == 0) {
            // 去除在中心的四面体，保留原本的中心四面体
            continue;
        }
        for(int j = 0; j < 4; j++) {
            tetrahedron_cells_expand.push_back(cell[j]);
        }
    }
}
